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ABSTRACT 

Aims. We study the properties of MHD turbulence driven by the magnetorotational instability (MRI) in accretion disks. To do this we 
perform a series of numerical simulations for which the resolution is gradually increased. 

Methods. We adopt the local shearing box model and focus on the special case for which the initial magnetic flux threading the disk 
vanishes. We employ the finite difference code ZEUS to evolve the ideal MHD equations. 

Results. Performing a set of numerical simulations in a fixed computational domain with increasing resolution, we demonstrate that 
turbulent activity decreases as resolution increases. The highest resolution considered is 256 grid cells per scale height. We quantify 
the turbulent activity by measuring the rate of angular momentum transport through evaluating the standard a parameter. We find 
a = 0.004 when (N,,Ny,N-) = (64, 100,64), a = 0.002 when (N,,Ny,N,) = (128,200, 128) and a = 0.001 when (N„Ny,N-) = 
(256, 400, 256). This steady decline is an indication that numerical dissipation, occurring at the grid scale is an important determinant 
of the saturated form of the MHD turbulence. Analysing the results in Fourier space, we demonstrate that this is due to the MRI 
forcing significant flow energy all the way down to the grid dissipation scale. We also use our results to study the properties of the 
numerical dissipation in ZEUS. Its amplitude is characterised by the magnitude of an effective magnetic Reynolds number Rbm which 
increases from 10"* to 10^ as the number of grid points is increased from 64 to 256 per scale height. 

Conclusions. The simulations we have carried out do not produce results that are independent of the numerical dissipation scale, 
even at the highest resolution studied. Thus it is important to use physical dissipation, both viscous and resistive, and to quantify 
contributions from numerical effects, when performing numerical simulations of MHD turbulence with zero net flux in accretion 
disks at the resolutions normally considered. 
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1. Introduction 

A long standing issue in accretion disk theory has been to iden- 
tify the source of anomalous transport of angular momentum. To 
date, the most likely mec hanism is believed to be t he magnetoro- 
tational instability (MRI: lBalbus & Hawlevli 19981) which simply 
requires a weak magnetic field and a radially decreasing angular 
velocity to operate in a highly conducting disk. Appropriate con- 
ditions are readily realised in many astrophysical accretion disks 
and the linear instability grows on dynamical timescales. Its non- 
linear evolution has been widely studied since the early 1990's. 
Local simulations using the sh earing box model dHawlev et al.l 
[199 5; Brandenburg^LalJll995|) were found to give rise to MHD 
turbulence with an associated rate of angular momentum trans- 
port compatible with the observ ations, with a, the standar d pa- 
rameter in standard disk theory (IShakura & Sunvaevlll973h . be- 
ing in the range 10^^-0.1 depending on the geometry of the mag- 
netic field. 

In this paper we focus on MHD simulations in a shearing box 
threaded by zero net magnetic flux initially and use the opera- 
tor split code ZEUS dHawlev & Stonelll995l). T he first detailed 
consideration of this case was bv lHawlev etd] (fT996). The po- 
tential importance for angular momentum transport in accretion 
disks is that it offers the possibility of local turbulence cou- 
pled with genuine dynamo action. If such activity can be main- 
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tained, it would be independent of any imposed magnetic field 
and being local, independent of distant boundary conditions. 
Accordingly this would be a robust outcome of the MRI, provid- 
ing a guaranteed level of tr ansport. However, some recent results 
dGardiner & Stonell2005b l) indicate that the saturated turbulent 
state is sensitive to numerical resolution. Accordingly, issues re- 
main as to whether a numerically converged saturated turbulent 
state can be achieved. It is important to note that should such 
simulations ultimately yield negligible or zero transport, the 
MRI can still operate to produce sustained turbulence and trans- 
port, but more attention would have to be paid to imposed fields 
and boundary conditions in the context of global simulations 
which have been and are being currently carried out CHawlevI 
l2001l:ISteinacker & Papaloizoul2002l; [Fromang & Nelsonl2006l) . 

The plan of the paper is as follows: In section[2]we describe 
the computational set up for the simulations we performed and 
algorithm used. We then go on to present results for typical runs. 
In section[3]we discuss the effect of resolution on the results and 
in section [4| we discuss the power spectra associated with the 
saturated turbulent states and use these to show that significant 
flow energy is always driven to the smallest numerically realis- 
able scales. Thus results remain dependent on resolution at the 
highest resolution studied. We then go on to discuss these results 
and their implications for understanding the non linear outcome 
of the MRI in section[5]and give our conclusions in section[6| 
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2. Initial conditions and run setup 

In this paper we solve the ideal M HD equations in a shearing box 
tooldrei ch & Lvnden-BeUlll965h using ZEUS ( iHawlev & Stond 
[1995). To do this we adopt a Cartesian coordinate system (x, y, z) 
with unit vectors (i, j, k) pointing in the directions of the coor- 
dinate axes. The vertical direction is defined by k and the az- 
imuthal direction by j. The frame rotates with the angular veloc- 
ity of a free particle in circular orbit at the centre of the box and 
the origin of the coordinate system £2 = Q.k. In this frame the 
basic equations can be written as 



dp 



0, 



dv 1 

+ p(v-V)v + 2p£l X V = -VP + — (VxB)xB, 

ot 4n 

dB 

— = Vx(vxB). 

ot 



(1) 
(2) 
(3) 



Here p stands for the gas density, v for the velocity, B for the 
magnetic field and P for the pressure. 

To close the system written above, one need to specify the 
pressure through an equation of state. For reasons of simplicity, 
throughout this paper we adopt an isothermal equation of state 
for which 



p-p4- 



(4) 



As usual, the ratio between the speed of sound cq and the angular 
frequency Q. can be used to define a disk scale height H. 

Given the above framework, a simulation is defined once 
the size of the box, the resolution and the initial geometry and 
strength of the magnetic field are chosen. In this paper, for the 
most part, we consider computational boxes of size (L^, L,., L,) = 
{H,nH,H), although we have also considered boxes of size 
(H, 2nH, H) in order to compare our simulations with those pub- 
lished in the literature previously. The resolutions used vary from 
(N^,Ny,N,) = (64, 100,64) to {N,-,Ny,N,) = (256,400,256). 
We will measure times in units of the orbital period, T - 2jt/Q.. 
As mentioned in the introduction, we focus exclusively on the 
special case in which no net magnetic flux threads the box (in 
the vertical or azimuthal directions) initially. The magnetic field 
at the start of the simulation is purely vertical and defined as 
follows: 



B- = Bosin(27rjc///) 



(5) 



where Bo is set such that the volume average ratio between ther- 
mal and magnetic pressure < P > equals 400. We checked, how- 
ever, that the saturated state of the turbulence depends neither 
on that value nor on the geometry of the field provided the net 
flux remains zero. At the beginning of each simulations, ran- 
dom velocity fluctuations of small amplitude are applied to an 
initial state with uniform gas density and velocity that is en- 
tirely due to the background Keplerian shear and takes the form 
V = (0, -3Q.X/2, 0). All the simulations performed here had the 

Courant number C - 1/2. 

Following standard practise (iHawlev et al.ll 19951) strictly pe- 
riodic boundary conditions are applied in y and z while boundary 
conditions that are periodic in shearing coordinates are applied 
in X. The latter require some care as they might introduce spu- 
rious numerical artefacts. Here, we applied these shearing box 
boundary conditions directly to the magnetic field. Although 
this procedure safely conserves the mean radial magnetic field 
threading the box to within round-off' error, the mean y and z 



components of the field are only conserved to within truncation 
error Thus there is a possibility of having long term accumula- 
tion of azimuthal or vertical magnetic fluxes, which would in- 
crease the MRI-induced turbulent activity in the box. In the fol- 
lowing, we will therefore monitor the time variation of the mean 
y and z components threading the computational domain during 
the simulations. 



2.1. Run parameters 

The details of the runs we performed are given in Table [T] The 
first column gives the simulation label. The second and third 
column give the box dimensions (Lj, L,., Lr) and the resolution 
(A^v, Ny, N^). The fourth column gives the simulation duration in 
orbital times. Finally, the last three columns give time averaged 
values of the Reynolds, Maxwell and total stresses, normalised 
by the initial thermal pressure Pq. They measure the rate of an- 
gular momentum transport and as usual are respectively defined 
through 



nRey 



1 



aRev 



CUM, 



= -^{P(VX - Vx)(Vy - Vv)) , 

^0 ^0 



J_l B,B, 
Po 



Po Po\ 47r 

a - QRey + auax - 



Po 



(6) 
(7) 
(8) 



In these equations, Vj and v,, are averages over y and z of the 
velocity components in the x and y directions respectively while 
the angled brackets denote a volume average. Note that because 
of the conservation of mass Po = (P). 

2.2. Standard runs 

To make a connection with previously published results, we first 
perform two runs with a moderate resolution. Both have 64 grid 
cells in the x and z direction. The first, labelled FS64, uses a box 
size (Lv, Ly, L^) - (H, 2nH, H) and 200 cells in th e y direction. It 
is alm ost identical to one of the runs presented bv lFleming et al.l 
(l2000h . The only difference is that these authors use an adiabatic 
equation of state and a somewhat lower resolution in the y di- 
rection. Model FS64 is compared with model STD64 in which 
the size of the computational box is halved in the y direction. To 
maintain the same effective resolution, A^y. is also decreased to 
100. Thanks to the improved computational resources that have 
become available in the last few years, models FS64 and STD64 
have been run for 300 and 1000 orbits respectively. 

The time history of aRey, auax and a are shown for both 
runs FS64 and STD64 in figure [T] respectively on the upper and 
lower panels. In each case, the dotted line represents URey, the 
dashed line shows umox while the solid line corresponds to a, 
the sum of the two. Both runs display the characteristic signa- 
tures of the MRI: an initial growth during the first few orbits due 
to the linear instability, a decrease of the stress after reaching a 
maximum as the linear instability breaks down into MHD tur- 
bulence and finally attainment of a saturated quasi steady state 
phase characterised by outward angular momentum transport for 
the remainder of the simulation. As observed in most simula- 
tions of this type, most of the transport is due to the contribu- 
tion of the Maxwell stress. Both runs show sustained MHD tur- 
bul ence for hundred s of orbits, in agreement with earlier stud- 
ies (ISano et alJl2004l) . with significant fluctuations occurring on 
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Table 1. Properties of the runs described in this paper: The first column gives the model label, while the next three columns give the 
size of the computational domain (Li, Ly,L ), the resolution (N_^,Ny,N:^) and the time (in orbits) for which the simulation was run 
respectively. The fifth to seventh columns indicate the rate of angular momentum transport by giving the volume and time averaged 
value of the Reynolds, Maxwell and total stresses, normalised by the initial thermal pressure, respectively (note that these values 
are obtained in each case by averaging the results from f = 40 orbits until the end of the simulation). 



0.020 



0.015 




time (orbits) 

Fig. 1. Time history of the stress parameters a/t^y, aMax and a 
for the runs FS64 {upper panel) and STD64 (lower panel). For 
each plot, the dotted curve corresponds to the Reynolds stress, 
the dashed curve corresponds to the Maxwell stress while the 
solid curve is the sum of the two. All of these are normalised by 
the initial thermal pressure. 



both short (less than 1 orbit) and long timescales (more than 10 
orbits). Time-averaged values of the stresses between f = 40 and 
the end of the simulation, are given for both runs in Table [T] For 
model FS64, we f i nd va lues almost identical to those reported by 
[Fleming & Stond (120031) in their Table 2. For model STD64, we 



find that the transport is weaker by about 30%, as for this model 
or = 4.1 X 10-^ whereas a = 5.9 x lO Mn model FS64. The dif- 
ference is due to the smaller box size in model STD64, as this is 
the only difference between the two simulations. Although this 
relation between the box size and turbulent activity is not yet 
well understood, it was already note d in earlier calculations of 
the shearing box jHawlev et alJll99ob . 

In conclusion, models FS64 and STD64 demonstrate good 
agreement between our results and previously published calcu- 
lations. The next step, which is the main goal of this paper, is 
to check the convergence of these results when resolution is in- 
creased. Given the large computational cost associated with well 
resolved simulations, it is necessary to choose one particular box 
size for these runs despite the differences between the runs FS64 
and STD64 mentioned above. To reduce the computational bur- 
den, we adopt the smaller box for the remainder of this paper, 
as the computing time is reduced by a factor of two. The con- 
sequences of changing the box size will be briefly considered in 
section |S!2l 



3. The effect of resolution 

In this section, we study the influence of resolution on the results 
of model STD64 by doubling (model STD128) and quadrupling 
(model STD256) the number of cells in each coordinate direc- 
tion. All other parameters are kept identical to those of STD64. 
The time history of the normalised stresses is shown in figure |2] 
for models STD128 and STD256 in the upper and lower pan- 
els (using the same conventions as in figure [T]i respectively. It is 
clear from figure|2]that aMax and ORe^ both decrease as resolution 
is increased. From Table[Tl the time averaged value of a in model 
STD128 is 2.2 x lO^^ while it is 1.1 x 10"^ in model STD256. In 
other words, in going from 64 grid points to 256 grids points per 
scale height, the turbulent activity decreases by approximately 
a factor of two each time the resolution increases by a factor 
of 2. Figure |2] also shows that the amplitude of the fluctuations 
of the stresses tends to decrease as the resolution is increased. 
This could be a signature of the decreasing importance of chan- 
nel flows that have been suggested as being responsible for these 
fluctuations (Sano''2007), as resolution is increased. 

Model STD256 is run for 105 orbits. The lower resolution 
models STD64 and STD128 show that this is enough to get a 
good estimate of the stresses. Indeed, averaging a in these two 
models between f = 40 and t - 105, we respectively found 
ff = 4.3 X 10-^ and 2.5 x 10-^ which are close to the values 
quoted in Table [T] and obtained by averaging over much longer 
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time (orbits) 



Fig. 2. Same as figure[T]but for the runs STD128 (upper panel) 
and STD256 (lower panel). Note the overall decreasing turbulent 
activity as the resolution is increased. 



periods. Averaging the stresses over 105 orbits should therefore 
be enough to get a good estimate of a. 

Next, we turn our attention to possible problems induced by 
the shearing box boundary conditions. To check whether y and 
z mean magnetic field are created in the box, we plot in figure [3] 
the time history of both in model STD128, normalised by Pq^- 
It shows no systematic increasing accumulation of net flux in the 
computational domain, despite the imperfect nature of the shear- 
ing box boundary conditions. Furthermore, the absolute value of 
both components is always very small. Their maximum strength 
during the simulation, expressed in terms of efifective beta values 
(defined as yS; - SttPq/ <Bi>^), is respectively = 5.4 x 10^ 
and ySj - 4.3 x 10'' for the y and z components. For the resolu- 
tion we are using, this is far too small a field strength to have 
any effect on the saturated state of the turbulence (the wave- 
length of the most unstable MRI mode for such weak fields is 
always smaller than a grid cell). We performed the same checks 
for models STD64 and STD256. For the former, the maximum 
values of the mean field components during the simulations cor- 
responded toPy = 2.3 X lO'' and ^, = 1.5 x 10^. For the later, 
we obtained Py = 3.8 x 10** and /?, = 1.7 x 10^ All these val- 
ues indicate that the boundary conditions have no effect on the 
results. 




100 150 
Time (in orbits) 



200 250 



Fig. 3. Time history of the mean azimuthal (solid line) and ver- 
tical (dashed line) magnetic field threading the box in model 

1 /2 

STD128, normalised by P^' . Both remain small enough no to 
aff'ect the long term evolution of the simulations. 




Fig. 5. Time history of the correlation length in the z direction 
of By for models STD64 (solid line), STD128 (dashed line) and 
STD256 (dotted line). As the resolution is increased, L,(By) de- 
creases, as might be expected from the appearance of the snap- 
shots shown in figure |4l 



Being confident that the models we present are not signif- 
icantly affected by the boundary conditions, we now turn to 
a more detailed analysis of their properties. To illustrate the 
changes in the structure of the flow as resolution is increased, 
figure |4] provides snapshots of the structure of the magnetic 
field. From left to right, contours of B,, in the (x,z) plane 
(y - 0) are given for models STD64 (left panel), STD128 (mid- 
dle panel) and STD256 (right panel). As the resolution is in- 
creased, smaller and smaller scale structure becomes apparent. 
The only limitation on the smallness of the scale appears to be 
due to finite resolution. It is possible to make this statement more 
quantitative by computing a vertical correlation length for By. 
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Fig. 4. Snapshots showing contours of the y-component of the magnetic field in the (x, z) plane (y - 0) for the runs STD64 (left 
panel), STD128 (middle panel) and STD256 (right panel). Smaller and smaller scale features in the magnetic field are seen as the 
resolution of the simulation increases. 



FollowinglL esur & Longarettil (l2007h . we define this correlation 
length through 



Lz(B,) 



i / / By(x,y = 0,z)By(x,y = 0,z')dz'dz 
[ jB^.(x,y^O,z)dz 



(9) 



In this definition, the symbol angled brackets denote an aver- 
age over the x direction. The time history of Ly(By) is plotted in 
figure |5] using a solid line for model STD64, a dashed line for 
model STD128 and a dotted line for model STD256. Because 
the raw data were very noisy, the initial curves were smoothed 
using a window of about 3 orbits, corresponding to 10 snap- 
shots. Figure |5] confirms the contraction of scale apparent from 
the trends seen in figure [U after about 40 orbits, L~(By) reaches 
a quasi steady value which decreases as resolution is increased. 
Averaging the results in time from f = 40 orbits until the end 
of the run, we obtained L^(By)IH = 0.06 for model STD64, 
0.04 for model STD128 and 0;025 for model STD256. For each 
model, these numbers indicate that typical size scale for struc- 
tures in By is a few grid cells. For model STD64, L,(By) corre- 
sponds to 3.8 grid cells, while it is equal to 5 grid cells for model 
STD128 and 6.5 grid cells for model STD256. Similar values are 
obtained when calculating a correlation length in the direction x, 
defined using an equation similar to Eq. (|9]l. 

All of these results indicate that the saturated state of MHD 
turbulence in these simulations is governed by the numerical dis- 
sipation of the code. It is therefore important to understand in a 
more detailed way the dissipation in ZEUS (and by extension 
any MHD code without specified diffusivities that relies on nu- 
merical dissipation to bound the size scale from below) and why 
and how it affects the results. In the following section, we use 
Fourier analysis to address this issue. 



4. Fourier analysis 

4.1. Power spectra 

Useful clues into the nature of MHD turbulence are usually pro- 
vided by power spectra. Here we compute the reduced power 
spectrum of kinetic energy in the vertical direction, which we 
define as 





E(k,) = ^po|v(fe,)|' 



(10) 



Fig. 6. Reduced power spectra of the kinetic energy (upper 
panel) and of the magnetic field (lower panel) in the z-direction. 
On both plots, the solid line corresponds to model STD64, the 
dashed line to model STD128 and the dotted line to model 
STD256. The dot-dashed line on the upper panel shows the slope 
A:""^^ expected in the standard Kolmogorov theory of incom- 
pressible hydrodynamic turbulence. In this and similar plots k is 
expressed in units of 1 /L. 



where \v(k,)\^ = |v,(/t,)P + \vy(k,)\^ + \v,(k,)\^. v~,(k,) is defined by 
vx(k,) jv,(x, y, z)e-"''~'dz >, (11) 
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where < . > stands for an average in the x and y directions. 
Similar definitions hold for v^ik,) and v,(L). In writing Eq. ( fTOl ). 
po stands for the (conserved) mean density of the flow. Similar 
expressions can be written to compute the reduced power spec- 
trum of magnetic energy. 

Both spectra are represented in figure |6] as a function of 
A:.. The upper panel shows the kinetic energy spectrum and the 
lower panel the magnetic energy spectrum. In both panels, the 
results of model STD64 are shown using the solid line, those of 
model STD128 are plotted with the dashed line and the dotted 
line finally represents the results of model STD256. At all reso- 
lutions, both kinetic and magnetic energy spectra show features 
typical of turbulence: the spectrum decreases with k,, showing 
that there is more energy at large scale. Note however the de- 
creasing power at large scales (both for the kinetic and magnetic 
energies) as resolution is increased. This is because turbulent ac- 
tivity (or, equivalently, angular momentum transport) decreases 
when resolution increases and is in agreement with the result of 
section[3] In the upper panel, the dot-dashed line enables the re- 
sults to be compared with the expected slope of a Kolmogorov 

— 11/3 

spectrum: E{k^ oc k, . There is no identifiable region with the 
expected Kolmogorov slope, that is maintained as the resolution 
is increased, that can be seen in the computed spectra, the best 
resolved calculation spanning almost two orders of magnitude in 
k,. For both the kinetic and magnetic energies, the spectra con- 
sist of a flat part at large scale, which grows in size as resolution 
is increased, and a decreasing part probably governed by numer- 
ical dissipation. In fact, these spectra fail to show any sign of an 
inertial range building up as resolution is increased. But is there 
any reason to expect these simulations to show a clear inertial 
range? Probably not. Because of the MRl, the flow is a priori 
unstable at all realisable scales and forcing and input is there- 
fore expected to occur all the way from the largest scale avail- 
able in the box down to the smallest MRl unstable scales (set 
by numerical dissipation). At these scales a small scale dynamo 
may al so operate and even transfer some energy back to larger 
scales (lBrummelletalJll998l: iBoldvrev et aLllioOSl: iPontv et alj 
l2005h . Thus in our case there is no good reason to suppose that 
any region of Fourier space is expected to be exclusively trans- 
ferring kinetic or magnetic energy downward to smaller scales, 
as would be required for an inertial range to be observed. To 
demonstrate this more clearly, we consider the properties of the 
Fourier transformed induction equation in the next section. 

In this context we comment that a situation where the MRl 
leads to dynamo activity differs from those such as occur when 
hydrodynamic phenomena such as the Rayleigh Taylor instabil- 
ity or Kelvin-Helmholtz instability produce turbulence. In the 
case of the Rayleigh Taylor instability, the source of energy is 
confined to the largest scales. Even though there are small scale 
instabilities in the linear regime, in the non linear regime these 
are overwhelmed by the advection process that results in the 
production of ev en smaller scales where dissipation takes place 
(IChertkovl2003h . In the case of the Kelvin-Helmholtz instability, 
the situation is similar with a source of instabilit y occuring only 
at large scales w hen there is no dynamo action dNepveul Il985t 
lRvuetal.ll2000l) . 



4.2. Equations 

We consider Eq. Q and decompose the velocity as the sum of 
the mean shear flow and the turbulent velocity field: 



where the mean shear flow is simply the y and z average of the y 
component of the velocity. 



Vshix) = Vshj 



(x, y, z)dydz 



(13) 



Using this decomposition, the right hand side of equation (O can 
be expanded and written as the sum of five terms: 

dB dB dV,h 

^ = -V./,— + B_,—^j-ivt ■ V)B - (V ■ v,)B + (B ■ V)v,(14) 

ot oy ox 

where the dependence of first two terms on velocity is through 
Vsh only. These describe advection by the mean flow and stretch- 
ing of the radial magnetic field lines by the background shear 
We take the Fourier transform of this equation and dot the re- 
sulting equation with the complex conjugate of the Fourier trans- 
form of B. Denoting the later by B*{k) and noting that 



B(Jt) = B{x)ex^(-ik-x)ct'x, 



(15) 



defines a finite Fourier transform, we obtain an equation govern- 
ing the evolution of the magnetic energy density in Fourier space 
in the form 



1 d\B(k)\^ 

2 di 
where 



- A + S + Thh + Tdivv + Th^, , 



A 
S 



-"Re 
+'Re 



B\k) 



Thb = -'Re 



-Re 



Th, = +Re 



dB :lr V ri 



ox 



[ivt ■ Y)B]e-'^-^d^x 



B^ik)- J J ^(V ■Vt)Be-'^ ^d^x 



B*{k) ■ 



/// 



[{B ■ V)v,]e-'^-^d^x 



(16) 

(17) 
(18) 
(19) 
(20) 
(21) 



where Re denotes the real part is to be taken. Note that we ap- 
plied the remap procedure described in iHawley et alJ (Il995h to 
account for the shear when computing the Fourier transform in 
the X direction. 

During the saturated phases of the previously described 
simulations, the time derivative of \B( k)\^ should vanish 
on average. In simulat ions of the MRl dHawlev et al.l 119951 : 
[Brandenburg et al]|1995 ) it is quite generally found that quan- 
tities are stretched out in the direction of the shear and thus the 
length scales in the x and z direction associated with the satu- 
rated state are significantly smaller than t he length scale in th e 
y direction. See for example Figure 4 of iHawlev et al.l (Il995h . 
Therefore we shall for simplicity consider only the particular 
Fourier plane defined by A:, = and then we get A - 0. Under 
these assumptions, Eq. ( fTSl l simplifies to 



S + Thb + Tdivv + Tbv - . 



(22) 



V = Vsh+ Vt 



(12) 



In this equation, S describes how the background shear creates 
the y component of the magnetic field, Thb is a term that accounts 
for magnetic energy transfer toward smaller scales, Tdiw is due to 
compressibility and Thv describes how magnetic field is created 
due to field line stretching by the turbulent flow. Each of these 
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Fig. 7. Plots of the functions S (solid line), Tti, {dashed line), 
TdivY {dotted line) and Tbv {dotted-dashed line) as functions of 
k for model STD256. These are averaged in time and Fourier 
space as described in the text. Note that S is positive at all scales, 
which is simply describing the production of By by the back- 
ground shear. 



Fig. 8. Same as figure |2l but for the poloidal part of the mag- 
netic energy. The term T,^, is positive at all scales, indicating 
that MHD turbulence is forced by the MRI from the largest scale 
available in the simulation box down to the dissipative scale. 



1,5x10"" 



terms now depends on the wavenumber k - {k^, k,). To improve 
the statistics, we average them on shell of given modulus k - \k\ 
as well as over time. 

Of course, in a real simulation of the type considered here, 
magnetic energy is damped by numerical dissipation. Therefore, f 
a more realistic equation than Eq. ( l22b would be I 

S +Tbt + Tdiy, + Tb, + £),„,„, = . (23) | 

where £),„,„, accounts for the numerical dissipation (of course 
'realistic' physical dissipation can be treated in a similar way, see 
paper II). In the following section, we study the balance between 
the various terms of Eq. ^ for models STD64, STD128 and 
STD256. 



4.3. Results 

As model STD256 is the most detailed of our runs in term of 
resolution, we start by describing the results we obtained in this 
case before comparing with the other simulations. As mentioned 
above, the results of model STD256 were averaged in time to im- 
prove the statistics. We used 20 dumps spanning about 60 orbits 
from t - 45 until the end of the simulation. Figure |7] plots the 
four terms appearing in Eq. (l22l i versus k. The solid, dashed, dot- 
ted and dotted-dashed lines respectively correspond to the terms 
S , Thb, Tdivv and Tbv Only the first is positive, while the other 
three terms are mainly negative, except for Tbb which is posi- 
tive at the largest scale of the box. The term Tdiw, accounting for 
compressibility, is also seen to reach significant values, proba- 
bly because of the presence a strongly nonlinear waves in these 
simulations (Gardiner & Stone 2005b; Papaloizou et al. 2004). 
The large and positive value of the S term simply shows that 
By is created at all scales by the background shear. For the MRI 
to be a proper dynamo, however, there has to be a way through 
which poloidal magnetic energy is regenerated from this toroidal 
field. To study that mechanism, we redo the analysis presented 
above but concentrate on the poloidal part of the magnetic field 
= {Bx, 0, B,) rather than on the full magnetic field B. In that 



Fig. 9. The solid line shows the variation of the numerical dis- 
sipation Z),^„„, = ~{Tbb ^divv "^bv^ "^i^h wavenumber k for 
model STD256. It is compared, using the dashed line, with the 
spectrum D^., that would result from a purely resistive dissipa- 
tion corresponding to a magnetic Reynolds number /Je^vf = 10^ 
(see the text for the definition of Rem)- The thin solid line sim- 
ply indicate the location of the zero on the plot. Note the good 
agreement between Z)^„„ and Dres at small scales, while they dis- 
agree significantly at large scales, showing that ZEUS numerical 
dissipation departs from a pure physical resistivity. 



case, both S and A vanish and under the assumption that MHD 
turbulence is in steady state, Eq. ( |23] l reduce to 



^ bh^ ^ divv 



bv num 



= 0, 



where we have now 



Tbb = 



B;{k) 



{V, ■ V)Bpe-'''-^dh 



(24) 



(25) 
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with corresponding expressions for the other terms in Eq ( l24l i. 
The results we obtained for model STD256 are plotted on fig- 
ure [8] with the same conventions as figure |7] They indicate 
that r^j, is positive for all k through turbulent velocity fluctua- 
tions creating poloidal magnetic field through field line stretch- 
ing, r^j, reaches its maximum at k ~ 40, which corresponds to 
roughly 1 /8th the size of the box (or 32 grid cells at this reso- 
lution). Nevertheless there are non negligible contributions from 
the largest scale available in the box all the way down tok ~ 200, 
which corresponds to only a few grid cells. This is an indica- 
tion that the MRI is forcing the flow at all available scales and 
explains why the power spectra presented in section 14.11 fail to 
show any inertial range. 

Another application of the above analysis is to provide in- 
formation about the dissipative properties of ZEUS. Indeed, be- 
cause of Eq. ( |23] l. the sum of the three terms represented on fig- 
ure [8] has to be balanced by the numerical dissipation. If that 
dissipation was exactly equivalent to a resistive process having a 
resistivity rj, its spectrum Dresik) would be given by 

D,,,(k) = r]k^\Bp(k)f . (26) 

In Figure|9]we plot Z)^,„„ {solid line) and -Ores (dashed line). In 
the later case we chose rj in order to get a good fit between both 
curves. At large k (or small scales), it is seen that -Dres is in 
good agreement with Z)^„,„. This fit can be used to estimate the 
numerical resistivity rj of the code (at this particular resolution 
for smaU scales and for this specific flow). This translates into an 
effective magnetic Reynolds number /Je^ = coH/rj which is of 
the order of 10^. This equivalence has to be used with caution, 
however, as both Z)^,„„ and -Df^s deviate significantly for val- 
ues of k smaller than 80. A contribution to this difference may 
arise from the poorer statistics available at the largest scale of 
the simulations (and particularly to non negligible contribution 
of the time derivative term in Eq. ( fT6] l at these scales), but it is 
also likely to be due to the fact that numerical dissipation can- 
not be simply described as arising from a diffusion process, at 
least at large scales. Also, it should be noted that the maximum 
amplitude of Z)^„„, (and therefore the location where most of the 
dissipation takes place) occurs atk ~ 70-100, at which point the 
term 7"^^, is still very significant (see figure [8]). This is why the 
saturated state of MRI driven MHD turbulence depends on res- 
olution. Based on this analysis, we would therefore predict that 
increasing the resolution by another factor of two would give a 
different saturated rate of angular momentum transport. 

In order to further investigate the effect of resolution, we re- 
peated the previous analysis for models STD64 and STD128. 
For model STD64, we used about 90 dumps regularly spaced in 
time during the 1000 orbits of the simulations to average the re- 
sults. For model STD128, we used 60 dumps that cover the last 
200 orbits of the simulations. The results are summarised in fig- 
ures[TO]and[TT|respectively. On both figures, the left panel is the 
equivalent of figure [8] These confirm the results obtained using 
model STD256. The MRI forces the flow at all available scales 
in the computational box. The right panels of figures [TO] and [TT] 
are equivalent to figure|9]but for models STD64 and STD128 re- 
spectively. Both confirm that small scale dissipation in ZEUS is 
similar to that provided by a physical resistivity, with magnetic 
Reynolds numbers of the order of 10"* and 3x10** respectively. 
However, as also seen for model STD256, numerical resistivity 
departs from a physical resistivity at large scale in a way that de- 
pends on resolution (note that the amplitude of the oscillations 
seen in plotting the numerical dissipation is reduced compared 
to model STD256, which is illustrating the fact that the statis- 
tics are improved when the simulation is integrated longer; as a 



1.2 




Fig. 12. Reduced power spectrum of the velocity and magnetic 
field averaged over time in model STD256. The plots represent 
the quantities A:?|B(fe,)|^ (solid line) and k^\v(k^)\^ (dashed line) 
which would measure the rate of dissipation for constant difFu- 
sivities. The dashed curve peaks at smaller A:, than the solid curve 
indicating that the hydrodynamic dissipation length is larger than 
its MHD counterpart. This in turn indicates that the numerical 
Prandtl number of ZEUS is larger than 1 . 




Fig. 13. Same as figure[T2] but for model STD64 (left panel) and 
STD128 (right panel). Both confirm that the numerical Prandtl 
number for ZEUS is larger than unity, in agreement with the 
results obtained with model STD256. 



result, the deviation of the numerical dissipation from a purely 
Laplacian process appears more solid). Both model STD64 and 
STD128 indicate that the numerical scheme in ZEUS is such that 
the numerical resisitivity could be negative at large scale. This is 
also suggested by model STD256 (see figure |9j although poor 
statistics make that conclusion less clear in that case. This pos- 
sible antidifFusive behavior of ZEUS was in fact already pointed 
out by Falle (2002) for ID shock calculations and has recently 
been observed to occ ur when studying the pro pagation of tor- 
sional Alven waves (iLesafire & BalbusI l2007h . It is not clear 
whether this is intrisic to ZEUS numerical scheme, to peculiari- 
ties introduced by the shearing box boundary conditions, or to a 
combination of both. 

5. Discussion 

5.1. The magnetic Prandtl number for ZEUS 

The previous analysis explains why numerical simulations per- 
formed with ZEUS fail to converge when the resolution is in- 
creased. It also leads to an estimate of the numerical resistiv- 
ity associated with such a calculation. However, the numerical 
effective kinematic viscosity v, associated with the momentum 



S.Fromang & J.Papaloizou: MHD turbulence in accretion disks. I. Convergence 



9 



5.0X10" - 




Fig. 10. The left panel is the same as figure [8] but computed using results of model STD64. Likewise, the right panel is the same as 
figure |9] applied to model STD64. The dashed line use a resistivity rj such that Rcm - lO'*. 




Fig. 11. Same as figure[TOl but for model STD128. On the right panel, the value of rj used to fit the numerical results is such that 
Rcm = 3 X lO"*. 



equation, has not been discussed. On general grounds one ex- 
pects this to be of the same order of magnitude as the numerical 
resistivity. To be more quantitative, it would be necessary to per- 
form an analysis similar to that presented in section l4~2l but for 
the evolution of the kinetic energy per unit mass, in order to esti- 
mate an effective v, which could in turn be used to obtain a mea- 
sure of the numerical magnetic Prandtl number Pm for ZEUS: 



Pm — — . 



(27) 



Such a procedure is however complex and beyond the scope of 
this paper. 

As an alternative, we obtain an indication of the value of Pm 
by considering = kl\B{k,J\^ and = kl\v(k^)\^. Here 
B{kz) is, to within an ignorable constant factor, the quantity as- 
sociated with the magnetic field corresponding to v(A;.) for the 
velocity field defined through equations ( fTOlfTTT i. It may also be 
found by evaluating the Fourier transform ( fTSI ) for kx - ky = 0. 

These quantities would be proportional to the numerical re- 
sistive and hydrodynamic dissipation for k^ - ky - 0, if the later 
could be described by a simple diffusion process. As we demon- 
strated in the previous section, this is not generally the case, but 
it appears to be reasonable for the smallest scales in the case of 



resistive dissipation. In order to obtain an indication of the effec- 
tive value of Pm, we shall make the very reasonable assumption 
that the numerical viscous dissipation at small scales can be de- 
scribed in the same way. This is expected because there are no 
strong shocks and the order of the finite difference scheme is the 
same for the induction equation and the equation of motion. Also 
there are no added hyperdiffusive terms, which would require a 
dependence on higher powers of k. in Dres and D,,,,. 

Both Dres and D,.,-, are plotted in figure [12] for model 
STD256 using solid and dashed lines respectively. Both curves 
are time averaged between f = 50 orbits until the end of the sim- 
ulation and are normalised by their maximum values. Figure [T2] 
shows that £),,„ peaks at a wavenumber k^^ which is smaller than 
the wavenumber k'l at which Dres reaches its maximum. These 
peaks should indicate the scale at which most of the dissipation 
occurs. In other words, they can be used for order of magnitude 
estimates of both the viscous and the resistive lengths ly and 
From figure[T2] we find fc'' ~ 100 and k'l ~ 150. This means that 
Pm is of order unity, and probably biased toward values larger 
than one. Indeed, since kl is smaller than k], it follows that the 
viscous length is larger than the resistive length, or that numeri- 
cal viscosity should be larger than numerical resistivity (see also 
the discussion in paper II). 
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me (orbits) 



Fig. 14. Time history of the stresses for the run STD64a. The dot- 
ted curve corresponds to the Reynolds stress, the dashed curve 
corresponds to the Maxwell stress and the solid curve gives the 
sum of the two. All of these are normalised by the initial thermal 
pressure. 




Time (orbits) 

Fig. 15. Time history of the Maxwell stress for the runs LB64 
(solid line) and LB128 (dashed line). For the former, time aver- 
age between t - 40 and the end of the run gives a = 5.2 x 10"^' 
while a = 3.2 x 10"^ for the latter Thus angular momentum 
transports decreases when resolution is increased. 



Of course, there is significant uncertainty associated with the 
above estimate and with the method we used to derive it, but 
the fact that Pm is larger than 1 appears to be solid. It is fur- 
ther confirmed by the results shown in figure [13] which is the 
same as figure[T2]but for model STD64 (left panel) and STD128 
(right panel). In both cases, the dashed curve peaks at smaller 
wavenumbers than the solid curve, indicating that the viscous 
length is larger than the resistive length, in agreement with the 
discussion above. 



5.2. Scaling arguments 

We note that the shearing box equations and boundary condi- 
tions may be transformed to a dimensionless representation. This 
is the case either, when the evolution is governed by partial dif- 



ferential equations, or by the finite difference equations of a nu- 
merical scheme. The transformation is performed by choosing 
L, the box size, Q ' and poL^ as the units of length, time and 
mass respectively. The resulting equations then depend only on 
the dimensionless quantities H/L, h/L and C, where h is the grid 
spacing and C denotes the Courant number (assuming a fixed 
aspect ratio for the grid cells). The unit of magnetic field is then 
'\JL^Q.^Po. Consequently we expect any one of the stress param- 
eters to have the scaling 



a oc (L/H)^ F(h/L, L/H, C), 



(28) 



where F is some unspecified function. For our calculations C - 
1/2 is fixed while runs STD64, STD128, and STD256 which 
have fixed L and H indicate that F is oc /i under those constraints 
and within the range of h considered. Then we may write 



a oc (Lh/H^)G(L/H). 



(29) 



For some unspecified function G. We have investigated the form 
of G by performing simulation STD64a. This has L reduced by a 
factor of two and h increased by a factor of two when compared 
to STD256. Thus if G(L/H) were constant, the stress parameters 
should be the same for the two runs. In fact the data show that on 
average the stresses were about 80 percent larger in STD64a (see 
figure[T4]and table[Tl respectively giving the stresses time histo- 
ries and averaged values). On the other hand the stresses showed 
significantly stronger time variability in that case but with a base 
level comparable to that in STD256. These results indicate that 
qualitative as well as quantitative changes occur when dimen- 
sionless parameters are varied. These differences may be due to, 
for example, a variation of the importance of compressibility as 
has been considered by Sano et al. (2004), or a variation of the 
effective Prandtl number The importance of the Prandtl number 
as determined by the physical diffusion coefficients when these 
determine the form of the saturated state is considered in a com- 
panion paper 

It is also of interest to ask whether the scaling of a with res- 
olution described in the present paper and expressed by Eq. ( |29] ) 
still holds for larger boxes. This could have important impli- 
cations for global simulations. Thus we performed two addi- 
tional simulations with a box of size (Lj , Ly, L,) - (2H, InH, H). 
The first one, labelled LB64, has a resolution (N,„Ny,N,) = 
(128, 200, 64), which amounts to 64 grid points per scale height. 
In the second, LB128, the resolution is doubled. FigurefTSjillus- 
trates the results through the time history of the Maxwell stress 
(the solid line corresponds to model LB64 and the dashed line 
to model LB 128). Again, we found a significant decrease of the 
turbulent activity as resolution is increased: a = 5.2 x lO""* for 
model LB64 and 2.8 x 10"^^ for model LB128. This tends to in- 
dicate that the results we present in this paper could extend to 
global disk simulations. Unfortunately the very high computa- 
tional demands associated with these simulations precludes ex- 
tensive studies at this time. 



6. Conclusion 

In this paper, we have shown that angular momentum transport 
induced by MHD turbulence decreases when the resolution is 
increased in numerical simulations performed with ZEUS in a 
shearing box in the absence of net magnetic flux. We have shown 
that this is due to the MRI forcing the flow at all scales, including 
those at which dissipation takes place. There is enough energy 
at these smallest scales to affect mean stresses in the saturated 
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state. All our results, taken together, demonstrate that it is impor- 
tant to use explicit diffusion coefficients that are large enough to 
produce more dissipation than numerical effects in local numeri- 
cal simulations of MHD turbulence with zero net flux performed 
with a finite difference code like ZEUS at cuiTently feasible res- 
olutions. Recent numerical simulations of MHD turbulence in 
the shearing box in the presence of an imposed magnetic flux 
showed that a also depend s on physical dissipation in that case 
dLesur & Longarettill2007l) . as it was found that it increases as 
the ratio of kinematic viscosity to magnetic diffusivity does. 

We note that there is no reason why this state of af- 
fairs shoul d not apply when usin g Godunov cod es like 
ATHENA jGardiner & Stond l2005al) or RAMSES ( Tevssieij 
I2OO2I; iFromang et al.1 1200617 or other codes making use of hy- 
perviscosity to stabilise the numerical scheme, like the PENCIL 
code ([Brandenburg & Dobler 2002) for example. 

A study of the effects of magnetic diffusivity and kinematic 
viscosity on local numerical simulations of MHD turbulence 
with zero net flux is the subject of a companion paper. 

Finally we comment that the results of this paper apply to the 
very simple computational set up of a local unstratified shearing 
box with zero net flux and for a restricted domain in parameter 
space. They should not be applied to more complex stratified or 
global simulations which will require separate studies. Neither 
should they be extrapolated beyond the parameter ranges con- 
sidered. 
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